function [LL, gradifo, grad] = projectAllAndLL( data, opts , prepgradhess)
% [LL, gradifo, grad] = projectAllAndLL( data, opts, prepgradhess )
% From all images (output of predictAllImagesfun)
% compute - log likelihood for each image.
% prepgradhess : switch that determines which info is put into gradifo.
%     0  : dont prepare gradient and/or hessian (into gradifo)
%     1  : just prepare gradient
%     2  : also prepare hessian.
%
% Created by Dirk Poot, Erasmus MC, 27-10-2011
if nargin<3 || isempty(prepgradhess)
    prepgradhess = 1;
end;
if nargout>=2
    prepgradhess = max(1,prepgradhess);
end;
szData = size(data);
if isempty(opts.project)
    if prepgradhess>0
        error('Apparently this case is (unexpectedly) used; but its not implemented yet. (I can only use global optimization in fit_MRI when project is used).');
    end;
    % make code correct:
    [lgpdf] = opts.logPDFfun( opts.data_in, data , opts.noiseLevel, [false true false]);
    LL = - sum( lgpdf(:) );
else
    datap = permute( data, [2:numel(szData) 1]);
    sel = repmat({':'},1,numel(szData));
    LL = zeros(size(opts.project,1),1);
    projgradifo   = cell( size(opts.project,1) , min(1,prepgradhess) );
    lgpdfgradhess = cell( size(opts.project,1) , prepgradhess );
    if nargout>=3
        grad = cell(1,size(opts.project,1));
    end;
    for k=1:size(opts.project,1)
        sel{end} = opts.projectSelectSource{k};
        if size(opts.noiseLevel,1)==1
            noiseLevel = opts.noiseLevel;
        else
            noiseLevel = opts.noiseLevel(k,:);
        end;
        if numel(noiseLevel)==1 && iscell(noiseLevel)
            noiseLevel = noiseLevel{1};
        end;
        [proj, projgradifo{k,:}] = project1(   datap(sel{:}) , opts.projectParameters{k} , opts.project(k,:), 1);
        if iscell( opts.logPDFfun ) || iscell(opts.project{k,1})
            % feed part of data and project to each logPDFfun.
            lgpdf = zeros(size(proj));
            lgpdfgradhessk = cell(numel(proj) , size(lgpdfgradhess,2));
            for pdfidx = 1:numel(lgpdf)
                if iscell( opts.logPDFfun )
                    pdffun = opts.logPDFfun{pdfidx};
                else
                    pdffun = opts.logPDFfun;
                end;
                if iscell(noiseLevel)
                    noiseLevelk = noiseLevel{1}(pdfidx);
                else
                    noiseLevelk = noiseLevel;
                end;
                [lgpdfk, lgpdfgradhessk{pdfidx,:}] = pdffun( opts.data_in{k}{ pdfidx }, proj{ pdfidx } , noiseLevelk, [false true false]);
                lgpdf(pdfidx) = sum(lgpdfk(:));
            end;
            for nargs = 1:size(lgpdfgradhess,2)
                lgpdfgradhess{k,nargs} = {lgpdfgradhessk{:,nargs}};
            end;
        else
            [lgpdf, lgpdfgradhess{k, :}] = opts.logPDFfun( opts.data_in{k}, proj , noiseLevel, [false true false]);
        end;
        LL(k) = - sum( lgpdf(:) );
        if nargout>=3
            % compute grad =  dLL/ d data_sel
            %          = dLL/d lgpdf * d lgpdf / dproj * dproj / d datap_sel
            % dLL/d lgpdf         = -1
            % d lgpdf / dproj     = lgpdfgradhess{1}
            % dproj / d datap_sel = project{k,3}( .., projgrad{1} )
            [grad{k} , projgradifo{k,1}] = project1( lgpdfgradhess{k,1} , projgradifo{k, 1} , opts.project(k,:), 3 );
            if numel(sel{end})==1
                grad{k} = reshape(-grad{k},[1 szData(2:end)]);
            else
                grad{k} = permute( reshape(-grad{k},[szData(2:end) numel(sel{end})]), [numel(szData) 1:numel(szData)-1]);
            end;
        end;
    end;
    if nargout>=2
        gradifo.projectGradinfo = projgradifo;
        %gradifo.lgpdfgrad = {lgpdfgradhess{:,1}};
        if prepgradhess>=2
            gradifo.lgpdfhess = {lgpdfgradhess{:,2}};
        end;
        if nargout>=3
            if isequal([opts.projectSelectSource{:}],1:szData(1))
                grad = vertcat(grad{:});
            else
                tmp = zeros(szData);
                for k=1:numel(grad)
                    tmp(opts.projectSelectSource{k},:) = tmp(opts.projectSelectSource{k},:) + grad{k}(:,:);
                end;
                grad = tmp;
            end;
        end;
    end;
end;